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ABSTRACT 

The transport of angular momentum by magnetic fields is a crucial physical process in formation 
and evolution of stars and disks. Because the ionization degree in star forming clouds is extremely low, 
non-ideal magnetohydrodynamic (MHD) effects such as ambipolar diffusion and Ohmic dissipation 
work strongly during protostellar collapse. These effects have significant impacts in the early phase 
of star formation as they redistribute magnetic flux and suppress angular momentum transport by 
magnetic fields. We perform three-dimensional nested-grid radiation magnetohydrodynamic (RMHD) 
simulations including Ohmic dissipation and ambipolar diffusion. Without these effects, magnetic 
fields transport angular momentum so efficiently that no rotationally supported disk is formed even 
after the second collapse. Ohmic dissipation works only in a relatively high density region within the 
first core and suppresses angular momentum transport, enabling formation of a very small rotationally 
supported disk after the second collapse. With both Ohmic dissipation and ambipolar diffusion, these 
effects work effectively in almost the entire region within the first core and significant magnetic flux 
loss occurs. As a result, a rotationally supported disk is formed even before a protostellar core forms. 

The size of the disk is still small, about 5 AU at the end of the first core phase, but this disk will grow 
later as gas accretion continues. Thus the non-ideal MHD effects can resolve the so-called magnetic 
braking catastrophe while maintaining the disk size small in the early phase, which is implied from 
recent interferometric observations. 

Subject headings: stars: formation - ISM: clouds - ISM: jets and outflows — radiative transfer — 
magnetohydrodynamics 


1. INTRODUCTION 

Circumstellar disks are supposed to form as natural by-products of star formation processes due to large angular 
momenta in natal molecular cloud cores. Although extensive observational and theoretical studies have been done, 
theories of disk formation are not complete yet. A disk plays crucial roles in star formation and evolution since its 
presence affects the evolution of the central star and its feedback. Moreover, circumstellar disks, or protoplanetary 
disks, are very important not only in the context of star formation but also of planet formation because they are the 
sites of planet formation. Therefore, it is of crucial importance to construct a consistent model of circumstellar disk 
formation in the context of star formation. 

The key to understand disk formation is angular momentum redistribution during protostellar collapse. Gravitational 
torque is important when magnetic fields are weak (Bate 1998; Saigo et al. 2008; Bate 2010; Tomida et al. 2010a; Bate 
2011), but typically magnetic fields observed in molecular clouds (Crutcher 2012; Li et al. 2014a, and references 
therein) are strong enough to remove angular momenta efficiently and overcome the centrifugal barrier (Mouschovias 
& Paleologou 1979, 1980; Basu & Mouschovias 1994, 1995a,b; Tomisaka 1998, 2000, 2002; Machida et al. 2005b,a; 
Commergon et al. 2010; Tomida et al. 2010b; Machida et al. 2011a; Bate et al. 2014). However, this magnetic braking 
is in fact too effective and strongly suppresses formation of rotationally supported disks and binary/multiples at least 
in the early phase of star formation and with the ideal MHD approximation, which is called the magnetic braking 
catastrophe (e.g. Mestel & Spitzer 1956; Mellon & Li 2008; Allen et al. 2003; Li et al. 2011, 2014b) or fragmentation 
crisis (Hennebclle & Teyssier 2008). 

Obviously this problem cannot be real because many protoplanetary disks and extra-solar planets have been observed, 
and because the fraction of binaries or multiples is known to be high (e.g. Raghavan et al. 2010; Duchene & Kraus 
2013). Therefore, there must be something missing in those early simulations which were highly idealized and simplified. 
This problem has been actively discussed recently and many solutions have been proposed. For example, non-ideal 
MHD effects such as Ohmic dissipation and ambipolar diffusion can redistribute magnetic flux and suppress angular 
momentum transport (Dapp & Basu 2010; Machida & Matsumoto 2011; Machida et al. 2011a; Dapp et al. 2012; 
Machida & Hosokawa 2013; Tomida et al. 2013). Mis-alignment between magnetic fields and the rotational axis can 
reduce the efficiency of magnetic braking and enable disk formation when magnetic fields are not too strong (Joos et al. 
2012; Krumholz et al. 2013). Or simply, a circumstellar disk can form later when the surrounding envelope almost 
disappears and magnetic braking becomes inefficient (Mellon & Li 2009; Machida et al. 2011b; Machida & Hosokawa 
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2013). Turbulence that ubiquitously exists in star forming clouds is another promising mechanism, but there are 
various interpretations how it circumvents the magnetic braking catastrophe; diffusion or reconnection enhanced by 
turbulence can redistribute magnetic fields (Lazarian & Vishniac 1999; Santos-Lima et al. 2012, 2013), incoherent 
magnetic fields are less efficient to transport angular momentum (Seifried et al. 2013, 2014), it can naturally induce 
mis-alignment between magnetic fields and rotation (Joos et al. 2013), or simply turbulence carries additional angular 
momentum locally. Given such many solutions proposed, the magnetic braking catastrophe is not so catastrophic as 
it sounds, but it still remains as quantitative questions; when and how is a disk formed, and how large (massive) is it? 

It is observationally well known that star forming clouds are strongly magnetized. Magnetic fields in dense molecular 
cloud cores are slightly weaker than or comparable to the critical field strength required to support those cloud cores 
(Crutcher 2012; Li et al. 2013). Observations of circumstellar disks around young stellar objects also have been actively 
carried out. Maury et al. (2010) performed interferometric observations toward Class-0 sources, and compared them 
with numerical simulations with and without magnetic fields. They did not detect largely-extended circumstellar 
disks, suggesting that the effects of magnetic fields are significant in formation and early evolution of circumstellar 
disks. Tobin et al. (2012, 2013) (see also Sakai et al. 2014) claimed that they found a large circumstellar disk with a 
radius of ~ 120 AU around a very young Class-0 object L1527 IRS, one of the best-studied examples. While this result 
was striking, Ohashi et al. (2014) observed the same object with higher resolution and sensitivity using the Atacama 
Large Millimeter/submillimeter Array (ALMA), and concluded that the disk radius should be about or smaller than 
60 AU. Clearly more detailed and extensive observations are needed, but overall, these observations imply that young 
circumstellar disks should be small, while they could form in the early phase of star formation and should grow later. 

Thus magnetic fields play crucial roles in star formation, and therefore their effects must be carefully clarified 
and quantified. Because the ionization degree in star forming clouds are extremely low, redistribution of magnetic 
Helds by non-ideal MHD effects is naturally expected and has a significant impact on angular momentum transport 
(Wardle & Koenigl 1993; Fiedler & Mousclrovias 1993; Basu & Mouscliovias 1994; Ciolek & Mouschovias 1994; Desch 
& Mouschovias 2001; Nakano et al. 2002; Tassis & Mouschovias 2007; Mellon & Li 2009; Kunz & Mouschovias 2009, 
2010; Li et al. 2011; Dapp et al. 2012). Non-ideal MHD effects are also important in the context of the magnetic 
flux problem, which means that typical magnetic flux in initial molecular cloud cores is far larger than that in formed 
stars by orders of magnitude and therefore substantial flux redistribution is required (e.g. Paleologou & Mouschovias 
1983; Braitlrwaite 2012). In this work, we investigate the non-ideal MHD effects and disk formation during protostellar 
collapse. For this purpose, we perform three-dimensional RMHD simulations of protostellar collapse including non¬ 
ideal MHD effects. In Tomida et al. (2013) (hereafter Paper I), we found that Ohmic dissipation can redistribute 
magnetic flux from the dense region and enable formation of rotationally supported disks around protostellar cores. 
In this paper, we extend our simulations to include ambipolar diffusion, which is more efficient in the lower density 
region. This paper is organized as follows. The equations and numerical methods are explained in Section 2, and 
the models used in the simulations are described in Section 3. The results of the numerical simulations are presented 
in Section 4. Section 5 is devoted to discussion and conclusions are summarized in Section 6 . In the Appendix, we 
describe the newly implemented solver for the non-ideal MHD effects and demonstrate its validity. 

2 . METHODS 
2.1. Basic Equations 

In order to study protostellar collapse, many physical processes must be properly taken into account, such as 
MHD, self-gravity, radiation transfer, chemistry and non-ideal MHD effects. In Paper I, we only considered Ohmic 
dissipation as non-ideal MHD effects. In this work, we extend our three-dimensional (3D) nested-grid simulations to 
include ambipolar diffusion, which is caused by decoupling between neutral and charged particles. We adopt the single 
fluid approach assuming strong coupling (Mac Low et al. 1995; Duffin & Pudritz 2008; Masson et al. 2012). Note 
that the Hall effect, another possibly important non-ideal MHD effect (Li et al. 2011; Braiding & Wardle 2012), is not 
included in this work. The governing equations we use are as follows: 
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where p, v, B,p, T g , e, E r , F r ,P r , $ are the gas density, gas velocity, magnetic flux density, gas pressure, gas temper¬ 
ature, total gas energy, radiation energy, radiation flux and radiation pressure tensor, respectively, c = 2.99792 x 
10 lo cms -1 is the speed of light, G = 6.673 x 10 -8 cm 3 g -1 s -2 is the gravitational constant, and a r = 7.5657 x 
10 -15 ergcm -3 K -4 is the radiation density constant. From top to bottom, these equations represent conservation of 
mass, the equation of motion, the induction equation including Ohmic dissipation and ambipolar diffusion, the gas 
energy equation, the solenoidal constraint, the Poisson’s equation of gravity, and the radiation transfer equations. We 
adopt the gray flux limited diffusion approximation for radiation transfer (Levermore & Pomraning 1981; Levermore 
1984). The tabulated Ohmic resistivity 770 and ambipolar diffusion rate 7 /a are calculated using a chemical network 
(see 2.3.2). For the Rosseland and Planck mean opacities, up and <jp, we combine three tables; Semenov et al. (2003), 
Seaton et al. (1994) and Ferguson et al. (2005). 


2.2. Overview of the Numerical Scheme 

In order to achieve the huge dynamic range to resolve star formation processes from the molecular cloud core scale 
down to the stellar core scale, we use the 3D self-similar nested-grid technique (Yorke & Kaisig 1995; Ziegler & Yorke 
1997; Machida et al. 2005b,a; Tornida et al. 2013). We solve these equations using the simple operator-splitting 
technique, using the same time steps over the whole nested-grid hierarchy (the so-called shared time stepping). Finer 
levels are generated at the center of the parent level so that the local Jeans length is resolved at least with 16 cells 
(Truelove et al. 1997; Commergon et al. 2008; Joos et al. 2012). We stop the simulations when this condition cannot be 
satisfied, as unphysical fragmentation may occur. Each level consists of 64 3 cells, and 14 levels are created by the end 
of the first core phase. Typical resolution in the first core (3-5 AU from the center of the cloud, which is covered with 
level 1=12) is Aa; ~ 0.14 AU or higher. The MHD part is solved using the HLLD approximate Riemann solver (Miyoslii 
& Kusano 2005) and the mixed divergence cleaning method (Dedner et al. 2002). The multigrid solver (Matsumoto 
& Hanawa 2003) is used for self-gravity. The time step for the MHD and radiation transfer parts is determined by 
the Courant-Friedrich-Levy (CFL) condition for the MHD part, and the radiation subsystem is solved implicitly. The 
implementation of newly-introduced non-ideal MHD effects is described in the Appendix. Detailed descriptions of 
other parts including microphysical processes are given in Paper I. 

Only in the model with ambipolar diffusion, we additionally introduce a density floor where magnetic fields are 
dominantly strong so that the time step does not become too small. When the gas density is so low that the plasma 
beta /3 = Pgas/Pm&g gets below 10 -3 , the gas density is increased so that /3 becomes 10 -3 while the gas temperature 
and velocity are unchanged. This floor works only in the very low density region in the infalling envelope just above 
the first core. We monitor the mass added by this density floor and confirm it is very small, about 5 x 10~ 5 M e in 
total, or less than 0.1% of the first core mass. Even with this density floor, the time step is limited by the Alfven 
speed in the low density region above the first core in the model with ambipolar diffusion, and as a result an excessive 
number of time steps are required. 

It is important to solve the energy equations including the effects of radiation transfer, although radiation does not 
drastically affect the evolution in the early phase of star formation compared to the previous simulations using the 
barotropic approximation (e.g. Bate 1998; Saigo et al. 2008; Machida et al. 2008a,b). Radiation heating and cooling 
are only crudely treated in the barotropic approximation, and heating by the non-ideal MHD effects is completely 
neglected. In reality, the thermal evolution varies depending on many factors, but the barotropic approximation 
simply assumes that the gas temperature is a function of the local gas density. Therefore, it cannot reproduce the 
realistic evolution and may affect the results artificially. For example, the stability of the first core depends on the gas 
temperature especially when it is gravitationally unstable (Stamatellos & Whitworth 2009). Also, the heating by the 
non-ideal MHD effects affects the structure of the first core (see Paper I), and the rates of the non-ideal MHD effects 
are sensitive to the gas temperature. Thus it is crucial to treat both dynamics and thermodynamics consistently. 


2.3. Microphysics 
2.3.1. Update on the EOS 


We update the equation of state from Paper I. Now the partition function for vibration of molecular hydrogen is 

1 


•^vib,H2 


1 - exp (%£) 
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where d v ib ~ 6330 K is the excitation temperature of vibration of molecular hydrogen 5 . This formula is better than 
the previous one used in Paper I because it does not diverge in the low temperature limit. However, this change has 


5 For example, NIST Computational Chemistry Comparison and Benchmark Database, NIST Standard Reference Database Number 
101, Release 16a, http://cccbdb.nist.gov/ 
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TABLE 1 

Summary of the initial model parameters and results 


Model 

f2 (xlO -14 s -1 ) 

B 0 (pG) 

Po 

OD 

AD 

TFc(yrs) 

Afpc 

Afp C /3>Fc(cgs) 

RSD Formation 

t?disk(AU) 

/ 

2.4 

20 

3.8 

N 

N 

770 

0.03 

2,500 

Not Formed 

- 

O 

2.4 

20 

3.8 

Y 

N 

1,050 

0.04 

7,200 

After Second Collapse 

1 

OA 

2.4 

20 

3.8 

Y 

Y 

2,750 

0.075 

37,000 

Before Second Collapse 

5 


Notes. The first five columns are the initial model parameters: the angular velocity, the magnetic field strength, the normalized mass-to- 
flux ratio, whether Ohmic dissipation and ambipolar diffusion are introduced. Other parameters are common: M = 1 M@, R ~ 8800 AU, 
p c = 1.2 X 10 -18 gem -3 , To = 10 K and A 2 = 0.1. The last five columns indicate the results of the simulations: the first core lifetime, the 
first core mass, the first core mass-to-flux ratio (not normalized), remarks on formation of rotationally supported disks, and the disk radius 
at the ends of the simulations. See the text for details. 

very minor impacts on the results of hydrodynamic simulations because contribution from the vibrational degrees of 
freedom becomes significant only in a high temperature, and most of hydrogen molecules are already dissociated there. 
Other parameters are unchanged including the ortho:para ratio of molecular hydrogen, which is 3:1. 

2.3.2. Ohmic Resistivity and Ambipolar Diffusion Rate 

The Ohmic resistivity and ambipolar diffusion rate are calculated based on the method described in Nakano et al. 
(2002) and Okuzumi (2009). This model solves a chemical reaction network of H 3 ", HCO + , He + , Mg + , C + , H + and 
dust grains, and derive equilibrium solutions. In addition, we consider thermal ionization of potassium (K), which 
produces a sudden increase of the ionization degree and sharp cut-off in the rates of the non-ideal MHD effects around 
T = 650 K. We adopt the reaction rates from the UMIST RATE 2012 database (McElroy et al. 2013). Dust grains can 
be charged from -14 to +14. We assume that the dust to gas ratio is 0.01, the dust density is 2gem -3 , and the size of 
the dust particles is 0.1 /.irn. The Ohmic resistivity in this work is higher than that in Paper I in the low density region, 
by about a factor of 3 in most regions except for the sharp increase around p ~ 10 - 14 gcm -3 (see Figure 3 below). In 
Paper I, this rise occurs around a higher density, p ~ 10 -13 gem -3 . These differences result in larger angular momenta 
remaining in the first cores, but the results are still qualitatively similar. 

We simply adopt a constant ionization rate 1.0 x 10 - 17 s -1 , assuming cosmic rays are the dominant ionization source 
and neglecting their shielding. In dense regions where cosmic rays cannot penetrate, this assumption overestimates 
the ionization rate and hence underestimates the diffusion rates. The ionization rate in such regions is generally 
determined by decay of radioactive elements. At least for our solar system, there is meteoritic evidence that short¬ 
lived radionuclides such as 26 Al existed when the solar nebula formed (Adams 2010). As discussed in Paper I, these 
short-lived nuclides yield an ionization rate of 0.7 — 1.0 x 10 - 18 s -1 (Umebayashi & Nakano 2009). In this case, 
neglecting cosmic-ray shielding would have only a moderate impact. By contrast, the ionization rate in dense regions 
can be extremely low if only long-lived radio active nuclides such as 40 K exist (Umebayashi & Nakano 2009; Kunz & 
Mouschovias 2009, 2010; Dapp et al. 2012). Because the lifetime of 26 Al is short and is comparable to the time scale 
of star formation, its abundance would vary in each system. Therefore this is an environmental parameter rather than 
an uncertainty. In this work, we aim to demonstrate that disk formation is possible even with the conservative case 
with the high ionization rate (i.e. the resulting non-ideal MHD effects are weak) as the first step. 

Another significant simplification in our model is that dust particles have a single population of 0.1 pm, neglecting 
their distribution and evolution. While dust grains can grow in collapsing clouds, they also can be spattered at shocks. 
On the other hand, dust grains evaporate at high temperatures. Although the effects of dust evaporation are taken into 
account in the opacities (Semenov et al. 2003), they are not included in our rates. Specifically, water ices evaporate 
around T ~ 150 K, while silicates evaporate around T ~ 1,400K. The latter has no significant impact on the rates 
because thermal ionization of potassium occurs below that temperature, and magnetic fields and gas are well coupled 
there. The evaporation of water ices, however, can be important but its effect strongly depends on dust properties 
such as composition and structures, which are highly uncertain. 

Thus the rates of the non-icleal MHD effects have large uncertainties, and they would vary in different environments. 
Therefore, despite the detailed calculation of the rates, our simulations should be considered as an experiment with a 
typical model, and broad parameter surveys must be performed in the future. 


3. MODELS 

We adopt the same initial conditions as in Paper I. An unstabilized Bonnor-Ebert (hereafter BE, Bonnor 1956; Ebert 
1955) sphere with m=2 perturbation is used as the initial density profile: 


p[r) 


PBE(r)(l + A 0 ) 


1 + A 2^ cos(2</>) 


(9) 


where pBE(r) is the the critical BE sphere with p c = 10 -18 g cm -3 and T = 10 K, R the radius of the critical BE sphere, 
A 0 the amplitude of the initial density enhancement, A 2 the amplitude of regularized to = 2 perturbation, respectively. 
We adopt Aq = 0.2 and A 2 = 0.1. The mass and radius of this cloud are M ~ 1M 0 and R ~ 4.25 x 10 -2 pc ~ 8800 AU. 
The initial free fall time at the center of the cloud is fg ~ 6.08 x 10 4 years. 

We initialize the cloud with solid-body rotation and uniform magnetic fields both aligned to 2 -axis. The angular 
rotation speed is D = 2.4 x 10 -14 s -1 or fits = 0.046. The magnetic field strength is B z = 20 p G and the corresponding 
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mass-to-flux ratio normalized by the critical value of stability for a uniform sphere is po — ^m/L) - t ~ 3.8 where 

$ = ttR 2 B 0 and (M/<I>) cr it = yjfr (^) (Mouschovias & Spitzer 1976). These parameters are the same to those in 
the fast rotating models RF ana IF in Paper I. 

The boundary conditions are also similar to those in Paper I. All the variables except magnetic fields outside the 
initial BE sphere are fixed to their initial values. Periodic boundary conditions are adopted for magnetic fields in order 
to avoid the divergence error. Because the gas outside the BE sphere is fixed, the magnetic fields are also anchored to 
the ambient gas. The boundary condition for self-gravity is set by the multi-pole expansion (Matsumoto & Hanawa 
2003). These model setups are chosen in order to model isolated dense molecular cloud cores. 

We calculate three models; Model I is the ideal MHD approximation, Model 0 includes Ohmic dissipation only, 
and Model OA includes both Ohmic dissipation and ambipolar diffusion. The model parameters are summarized in 
Table 1. While Model I evolves without any trouble until the central temperature exceeds T c ~ 80,000 K, the other- 
two models are more troublesome. As a result of suppression of angular momentum transport, rotationally supported 
disks are formed in these two models. These disks become gravitationally unstable after the second collapse in Model 
O and during the second collapse in Model OA. Because of the limitation of the nested-grid code, we only can refine 
the center of the computational domain, and therefore we cannot follow further evolution when fragmentation occurs 
outside the finest grid. Thus we stop the calculations around T c ~ 30, 000 K in Model O and T c ~ 2,100K in OA. 
For this reason, we focus on the first core phase in this work. Larger grids or more flexible adaptive mesh refinement 
(AMR) is required for further studies after fragmentation occurs. 

4. RESULTS 

4.1. Overview of Evolution 

Here we summarize the evolutions of the three models. The time evolution of the density and temperature at the 
center of the clouds as well as the first core masses are shown in Figure 1. Figure 2 is the evolution tracks of the 
central gas element in the p-T plane, and Figure 3 shows the rates of Ohmic dissipation and ambipolar diffusion at the 
center of the cloud in Model OA. The profiles of the density, temperature, velocities and various forces at the end of 
the first core phase (we simply define it as T c = 2, 000 K although the effective adiabatic index gets below the critical 
value 7 = 4/3 slightly below this temperature) are shown in Figure 4. Model I is essentially the same as Model IF in 
Paper I, and Model O is similar to RF in Paper I. 

Initially the clouds collapse isothermally because radiation cooling is very effective, but the gas starts to evolve 
almost adiabatically when the central density exceeds about p c ~ 10 -13 gcm -3 . Then the temperature quickly rises 
and quasi-hydrostatic objects, the first cores, are formed (Larson 1969; Masunaga & Inutsuka 2000). The formation 
of the first core is slightly earlier in Model OA, because ambipolar diffusion redistributes the magnetic flux before the 
gas becomes dense enough to form a first core, and it slightly weaken deceleration of the infall motion by magnetic 
fields. However, this is not significant because the dynamical collapse is much faster than the redistribution of the 
magnetic fields, and the difference is only 1% compared to the time from the beginning of the collapse to the formation 
of the first core. It should be noted that this difference depends on the stability of the initial condition. If the initial 
condition is sufficiently unstable (super critical), ambipolar diffusion has only minor impact on the initial isothermal 
collapse (see also Chen & Ostriker 2014; Heitsch & Hartmann 2014). 

Substantial differences arise during the first core phase. The first core evolves more slowly in the non-ideal MHD 
models, especially with ambipolar diffusion. The non-ideal MHD effects redistribute magnetic fluxes from the first 
cores and suppresses magnetic angular momentum transport. The additional rotational support changes the structures 
of the first cores and extends their lifetimes. Here we define the first core lifetime simply as the time between the epochs 
when the central temperature exceeds 100K and 2,000K. The lifetimes of the first cores are about 770 years in Model 
I, 1,050 years in O, but 2,750 years in OA. The first core in OA is qualitatively different and lives significantly longer 
than the others. While the accretion rates are essentially the same in all the models (initially M~4x 10 -5 M 0 yr _1 
and gradually decrease later), the final masses are significantly different, about 0.03M 0 in Model I, 0.04M 0 in O, 
and 0.075 M 0 in OA. The first core in the ideal MHD model is essentially the same with spherical first cores without 
rotation, while the first core in Model O is slowly rotating but still looks almost spherical as rotational support is 
not dominant (Figure 4, see also Paper I). The first core in Model OA becomes disk-like and is clearly supported 
by rotation, although contribution of the gas pressure is still substantial. The non-ideal MHD effects also provide 
additional heating as in eq.(4), and the gas element follows a hotter (i.e. higher entropy) evolution track in Model O, 
and it is more prominent in OA (Figure 2). These variations in the thermodynamic evolutions cannot be reproduced 
with the barotropic approximation. Because the hotter gas has higher pressure and can support more mass, this also 
contributes to the larger masses and longer lifetimes in the non-ideal MHD models. 

When the temperature reaches about 1,400K, all the dust components evaporate and the opacities drop drastically 
(Figure 2, see also Paper I). This effect is not obviously visible in the evolution track of Model / because it collapses 
very quickly. Model O, on the other hand, experiences considerable radiative cooling after dust evaporation and its 
evolution track merges to that of the ideal MHD model. In Model OA, the dust evaporation affects the thermal 
evolution more prominently because it evolves more slowly and therefore has more time to cool (see also Tonfida 
2014). It evolves along the dust evaporation temperature for a while 6 , and the gas gets colder than those in Model / 
and O. 

6 Although this behavior is interesting, this might be artificial because we use the tabulated dust opacities as functions of the density 
and temperature, assuming there is no hysteresis. Specifically, the evaporated dust grains immediately form again once the gas cools down 
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Fig. 1. - Time evolution of the gas density (top) and temperature (middle) at the center of the clouds, and the first core masses (bottom). 
The temperatures where the dust evaporation completes and where the second collapse begins are also shown with the gray dotted and 
dash-dotted lines. Note that these temperatures depend on the gas density, but the dependency is weak and the densities where these 
events happen in these models are close enough (see Figure 2). 

Dissociation of molecular hydrogen becomes significant when the temperature at the center of the clouds reaches 
about T c ~ 1, 800 K. Because of this strongly endothermic reaction, the gas pressure fails to maintain balance with the 
gravity and starts to collapse again. This second collapse is fairly violent in Model I and O. They dynamically collapses 
almost at the time scale of free-fall. After the second collapse completes, another hydrostatic objects, protostellar cores 
are formed. While the protostellar core in Model I is almost spherically symmetric, that in Model O is rotating and 
a disk is simultaneously formed (see also Paper I). Because the resistivity in this work is higher than that in Paper 
I, the disk is larger but its radius is still small at the end of the simulation, which is only R ~ 1 AU. On the other 
hand, Model OA evolves again differently; because of the centrifugal force, the second collapse occurs rather gradually. 
Radiation cooling is still effective even in this phase and the gas evolution follows the colder track, deviating from 
the adiabatic evolution (Figure 2). As mentioned above, we stop the calculations of Model O and OA earlier because 
fragmentation occurred. 

The evolution of the diffusion rates and corresponding magnetic Reynolds numbers evaluated at the center of the 
cloud in Model OA is shown in Figure 3. The Ohmic dissipation rate in Model O behaves almost the same. The 
ambipolar diffusion rate is higher than the Ohmic dissipation rate in almost all the regions, and Ohmic dissipation 

below the evaporation temperature. In reality, it would take some time to reform dust grains and their structure and/or size distribution 
would be different from those before the evaporation. While this reformation of dust grains is highly uncertain, fortunately, it would not 
affect the evolution significantly. The gas evolves along the dust evaporation line if the opacities after reformation are higher. On the other 
hand, it evolves almost isothermally until the gas opacities and density become large enough if the opacities are lower, but the evolution is 
already close enough to isothermal. 
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Fig. 2.— Evolution tracks of the gas elements at the center of the clouds in the p-T plane. The small inset shows the region near the 
onset of the second collapse. The gray-dotted line indicates the temperature at which dust grains evaporate completely (Semenov et al. 
2003). The evaporation begins below this temperature and the opacity is already reduced. Also, the dash-dotted line roughly shows the 
temperature where the second collapse begins (i.e. 7 4/3, see also Tomida (2014)). 



Fig. 3. — The rates and corresponding Reynolds numbers of Ohmic dissipation and ambipolar diffusion at the center of the cloud in 
Model OA as functions of the gas density. The magnetic Reynolds numbers are calculated assuming that typical velocity and length scale 
are the free-fall velocity and Jeans length; Rq a = r’frA..I /rjo A (see Paper I, Appendix C). 


becomes dominant over ambipolar diffusion only in the dense region where p c A) 2 x 10 10 gem 3 . The Ohmic dissi¬ 
pation rate is very low in the low density region, and it becomes active only in the dense region within the first core 
(2x 10 -11 genr -3 p c 6x 10 -10 gem -3 ). The ambipolar diffusion rate is much higher than the Ohmic dissipation rate 
in the low density region, but it does not change the dynamics substantially until the first core is formed. Ambipolar 
diffusion works effectively in the broader range than the Ohmic dissipation (5 x 10 -12 g cm -3 p c 5 x 10 -10 gcm -3 ), 
removes magnetic flux and suppresses angular momentum transport more significantly. Both the Ohmic dissipation 
and ambipolar diffusion rates drop sharply when the gas temperature exceeds about 650 K (or p c ~ 5 x 10 -lo gcm -3 ) 
because the gas and magnetic fields recouple by thermal ionization of potassium (K). 

4.2. Outflows 

First, we compare the outflows driven around the first cores. Vertical cross sections of the density and temperature 
in this scale are presented in Figure 5. These outflows are driven by the magneto-centrifugal force (Blandford & 
Payne 1982) and remove angular momenta from the central objects efficiently (e.g. Tomisaka 1998, 2002; Tomida et al. 
2013). The driving regions are located around R = 10 — 30 AU, corresponding to the outer-most peaks of the rotational 
velocity v^, in Figure 4. The gas density at these radii is still low and neither Ohmic dissipation nor ambipolar diffusion 
works strongly yet. The difference of the driving radii between the models is simply originated from the ages of the 
systems, because accreting angular momenta get larger as the systems evolve. 
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Radius [AU] Radius [AU] Radius [AU] 

Fig. 4.— Radial profiles of the gas density, temperature along the x- (in the disk mid-plane, red) and 2 -axes (along the rotational axis, 
blue), the infall (red) and rotation (blue) velocities along the ai-axis (in the mid-plane), and the radial accelerations by the various forces 
(The gas pressure gradient (red), the magnetic pressure gradient (green), the magnetic tension (blue), the centrifugal force (magenta), and 
the gravity (black). The radiation force is not plotted because it is negligibly small here) in the mid-plane at the end of the first core phase, 
from top to bottom. From left to right, the different columns are for Model 7, O and 0/1, respectively. 


The sizes of the outflows are significantly different between the models. Before the end of the first core phase, the 
outflows extend vertically about 70 AU in Model I, 80 AU in O, and 170 AU in OA. However, the structures and 
velocities of the outflows are similar, about 1 kms -1 in all the cases. This is because the non-ideal MHD effects, both 
Ohmic dissipation and ambipolar diffusion, only work in the high density region within the first core, and do not affect 
the driving regions which are located outside the first cores as discussed above (Figure 4). Therefore, the larger sizes 
in the non-ideal MHD models are solely consequences of the longer first core lifetimes by rotational support. 

While the outflow structures are quite similar, the temperature distribution in this scale is prominently warmer in 
Model OA than the other models. This is because the first core in OA is more massive and hotter, and heats up the 
surrounding material by radiation more strongly. Because of the longer lifetime and higher luminosity, the probability 
of finding a first core by observations would be higher than previously expected, although it should be still very rare 
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Fig. 5.— Vertical cross sections of the gas density (top) and temperature (bottom) in the outflow scales (1 = 8 or ~ 140AU for Model I, 
and l = 7 or ~ 280AU for the rest) at the end of the first core phase (T c = 2, 000 K). Projected velocity vectors are overplotted with white 
arrows. Because the outflows are roughly axisymmetric in the large scale, only vertical cross sections are shown here. 

because the lifetime is still much shorter than the timescale of the whole star formation process. 

Another noticeable difference in the outflows is that a polar cavity with a very low density is formed in Model OA. 
The gas near the pole is removed by accretion and the centrifugal force, rather than expelled by the outflow itself. The 
opening angle of this cavity is still very narrow in the early phase, but it will expand as the gas in the outer envelope 
with higher angular momenta will accrete later. Although it is not very likely, this cavity might enable us to observe 
the central object more directly if it is in a face-on configuration. 


4.3. First Cores and Second Collapse 

As we can see in Figure 3, the non-ideal MHD effects are most active in the dense region within the first cores and 
significant magnetic flux loss occurs where the Reynolds number gets below unity. As a result, much more prominent 
differences between the models arise in the first core evolution and structures. Figure 6 shows evolution of the masses, 
angular momenta and mass-to-flux ratios of the first cores. In order to produce these graphs, we define the first core 
by a simple criterion; where the gas density is higher than a critical density 7 p C rit = 10~ 12 gcnW 3 . Also, the density 
and temperature cross sections of each model are shown in Figures 7-9. Three epochs are shown in these Figures; 
just before thermal ionization of potassium (T c = 600 K, left), beginning of dust evaporation (T c = 1,400K, middle) 
and the end of the first core phase, just after the onset of the second collapse (T c = 2,000K, right). The magnetic 
field lines at the end of the first core phase are presented in Figure 10. The properties of the first cores and disks are 
summarized in Table 1. 


4.3.1. Ideal MHD Model 

The first core in Model I is essentially the same as that in Model IF of Paper I. Magnetic braking is so efficient that 
almost all the angular momentum in the first core is removed, and the first core is virtually spherically symmetric and 
not rotating (Figures 4 and 7). The first core mass and radius at the end of the first core phase are about O.O3M 0 
and 3 AU. The first core and outer pseudo-disk exhibit strong warping, which is driven by the magnetic interchange 
instability (Spruit et al. 1995; Stehle & Spruit 2001; Krasnopolsky et al. 2012). The angular momentum in the first core 
decreases as it evolves. The two sharp “wedges” are induced by strong perturbation from the magnetic interchange 
instability, indicating that the first core angular momentum is very small (essentially zero) and is easily dominated 
by the external warping. Because the thermal pressure and rotational support are smaller in this model than in the 
others, and the magnetic pressure is not significant, the first core mass is the smallest in this model. 

7 This value is one order of magnitude higher than the critical density used in Paper I. We have changed the threshold so as to exclude 
the high density region in the pseudo-disk (a flattened disk-like structure but not supported by gas pressure or rotation) especially in 
Model OA. This is why the angular momentum of Model I shown in Figure 6 behaves differently from that in Paper I, but in fact they are 
consistent if the same threshold is adopted. 
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Gas Density [g cm" 3 ] 


Fig. 6. — Evolution of the masses, angular momenta, and mass-to-flux ratios of the first cores from top to bottom as functions of the gas 
density at the center of the cloud. The red dotted line in the bottom panel is extrapolation of the mass-to-flux ratio before the magnetic 
interchange instability develops (see the text). Note that this mass-to-flux ratio is presented in the cgs unit and not normalized. 


To see the degree of magnetization, we look into the mass-to-flux ratio Mp c/^fc of the first cores. The magnetic 
flux threading a first core is measured by integrating the vertical magnetic flux density over the mid-plane as follows: 

$fc= [ B z dS= B * AS ’ (10) 

•tp>Pcrit p>Pcrit 

where AS = AxAy is the surface area of a discretized cell, and the summation is taken over the mid-plane. This 
formula captures the total magnetic flux correctly when the mid-plane is the plane of symmetry. In Model I , the 
symmetry breaks down by the magnetic interchange instability after the central density exceeds a few xl0 _9 gcm~ 3 , 
while the other two models maintain the symmetry during the first core phase. Therefore, the mass-to-flux ratio in 
Model I after the sharp increase is not accurate, and it only gives the upper limit. Instead, it is more appropriate to 
assume this value is almost unchanged beyond this point (the red-dotted line in Figure 6) because the time after this 
is very short and the first core acquires only a little mass (see Figure 1). Taking this point into account, the first core 
magnetization remains almost constant during the first core phase in the ideal MHD model. 

In the large scale, the magnetic fields look like the so-called hourglass shape. The field lines are tangled near and 
within the first core because of the magnetic interchange instability (Figure 10). Because the first core is not rotating, 
the field lines are not wound up. 

After the second collapse begins, the gas collapses dynamically and a protostellar core is formed very quickly in a few 
years. The protostellar core is also virtually spherical because almost no angular momentum exists in the protostellar 
core (at least shortly after its formation). The protostellar core properties are essentially the same as a spherical model 
without rotation and magnetic fields (Paper I). 
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Fig. 7.— Vertical (top two rows) and horizontal (bottom two rows) cross sections of the gas density (the first and third rows) and 
temperature (the second and fourth rows) in the first core scale (1=11 or ~ 18AU) of Model /, when the temperature at the center of the 
cloud is 600K (left column), 1,400K (middle) and 2,000K (right), tpc is the a g e of the first core after its formation (when T c = 100 K). 


log (Gas Temperature [K]) log (Gas Density [g cm' 3 ]) log (Gas Temperature [K]) log (Gas Density [g cm -3 ]) 
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4.3.2. Non-Ideal MHD Model with Ohmic Dissipation Only 

Model O is very similar to Model RF in Paper I, although the Ohmic dissipation rate in this work is higher than 
the one used in Paper I. Distribution of dissipative regions are shown in Figure 11 (the top row) using the magnetic 
Reynolds number Rq = |v|Aj/?7o where Aj is the local Jeans length. Because Ohmic dissipation works only in the 
high density region within the first core, the first core properties (Figure 6) are identical to those in Model I when 
it is formed. As the first core evolves, the magnetic flux is redistributed outwardly by Ohmic dissipation. The gas 
and magnetic fields recouple in the central region after the thermal ionization of potassium begins, but the outer 
region remains resistive and the mass-to-flux ratio increases monotonically. At the end of the first core phase, the 
mass-to-flux ratio is about three times higher than in Model I. Angular momentum transport by magnetic fields is 
significantly suppressed, and the first core acquires a large angular momentum, which also monotonically increases. 
As a result, the first core has considerable rotation and evolves slightly more slowly than Model I (Figure 1). However, 
the rotation is not fast enough to support the first core and its structure remains almost spherical (Figures 4 and 8). 
At the end of the first core phase, the mass and radius are about 0.04 M 0 and 5 AU, which are larger than the ideal 
MHD case because of additional support by rotation and heat produced by Ohmic dissipation (Figure 2, see also Paper 
I). However, these differences in the first cores are rather minor. Thus, roughly speaking, Ohmic dissipation does not 
have a drastic impact on the first core structure and evolution. 

The magnetic fields in the outer region of the first core are straightened out by the strong Ohmic resistivity. On 
the other hand, the magnetic fields are tightly wound up and amplified by the shearing rotation in the hot recoupled 
region near the center of the first core (Figures 10 and 11). The magnetic pressure in this region is not strong enough 
to launch outflows during the first core phase. Later in the protostellar core phase, strong toroidal magnetic fields are 
produced by the faster rotation in the disk and drive fast jets by the magnetic pressure (see Paper I for the detailed 
discussion about the protostellar cores). 

The second collapse occurs dynamically at the beginning as in Model I because the rotational support is not 
significant. However, since the gas spins up as it collapses, a rotationally supported disk forms around the protostellar 
core essentially in parallel. The disk size is very small, only about 1 AU at the end of the simulation (when T c = 
30,000 K). At the same time, fast bipolar jets are launched by the magnetic pressure of the toroidal magnetic fields 
amplified by the shearing rotation (Paper I, see also Machida et al. 2008a). This disk is larger than that in Paper I 
as a consequence of the higher Ohmic resistivity, but qualitatively the same. In this sense, Ohmic dissipation resolves 
the magnetic braking catastrophe and enables formation of a rotationally supported disk in the early phase of star 
formation. We stop the simulation because the Jeans condition cannot be satisfied beyond this point due to the 
limitation of the simulation code. 

4.3.3. Non-Ideal MHD Model with Ohmic Dissipation and Ambipolar Diffusion 

The first core in Model OA with both Ohmic dissipation and ambipolar diffusion is qualitatively different from 
the other two cases. Ambipolar diffusion is more effective than Ohmic dissipation especially in the low density gas 
(Figure 3), and suppresses angular momentum transport more significantly. The first core is already less magnetized 
from the beginning, and it loses the magnetic flux continuously during the first core phase (Figure 6). Figure 11 
indicates that the non-ideal MHD effects work strongly in almost the whole first core, especially in the early phase 
(see also Susa et al. 2015). At the end of the first core phase, the mass-to-flux ratio is higher than the ideal MHD 
model by more than one order of magnitude. 

Angular momentum transport by magnetic fields is strongly suppressed in such a diffusive region, and substantial 
angular momentum remains within the first core. As a result, the first core becomes a disk supported by the centrifugal 
force, although the gas pressure is still substantially contributing (Figure 4). The rotational profile is not Keplerian 
yet (i.e. the rotation profile is not v</, oc A -1 / 2 ) because of this gas pressure support and the self-gravity. The disk 
radius is initially about 6 AU when it is formed, gradually shrinks to ~ 5AU, and then remains almost constant, 
because magnetic braking is still working outside the resistive region and regulates angular momentum accretion. 
Although the disk remains small until the end of the simulation, it accumulates large angular momentum as the gas 
accretes from the envelope. Eventually it acquires almost one order of magnitude larger angular momentum than 
that in Model 0, while the mass (0.075M e at the end of the first core) is larger only by a factor of two. This fast 
rotation is sufficient to the trigger the gravito-rotational instability (Toornre 1964; Bate 1998; Saigo & Tomisaka 2006; 
Saigo et al. 2008; Tsukamoto & Machida 2011), and the disk becomes non-axisynnnetric in the late phase (Figure 9). 
Angular momentum transport by gravitational torque via non-axisynnnetric structures is usually less effective than 
magnetic braking, but it becomes important where the magnetic fields are weakened by the non-ideal MHD effects. It 
is naturally expected that this disk will grow later and the rotational support will be more significant as the gas with 
higher angular momenta will accrete from the outer region. 

The magnetic flux expelled by the non-ideal MHD effects piles up outside the diffusive region (R ~ 5AU), i.e. 
the first core. In this situation, the so-called magnetic wall can be formed, where the magnetic pressure strongly 
decelerates the infalling gas and possibly form a shock (e.g. Li & McKee 1996; Tassis & Mouschovias 2005a,b, 2007; 
Kunz & Mouschovias 2010). However, such a structure is not prominent in this model, although the accretion velocity 
is indeed slower than the free-fall velocity due to the magnetic fields (Figure 4). The outermost deceleration (at 
R ~ 100 AU) is formed by interacting with the bipolar outflow (see Figure 5), the second one (at R ~ 30 AU) is a 
centrifugal barrier, and the innermost one (at R ~ 5 AU) is the first core shock. This is likely because our rates of 
the non-ideal MHD effects are not high enough, or because the simulation duration is not long enough to accumulate 
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T c =600 K (t FC =570 yrs) 


T c =1,400 K (t FC =900 yrs) 


T c =2,000 K (t FC =1,050 yrs) 




Fig. 8 .— Same as Figure 7 but of Model O. 
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Fig. 9.— Same as Figure 7 but of Model OA. 
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sufficient magnetic flux. 

As discussed before, this model suggests that first cores are more readily observed because of the higher luminosity 
and longer lifetime. Another interesting observational aspect is its variability. While mass and angular momentum 
transport by magnetic fields are almost steady and smooth, the gravitational torque via non-axisymmetric structures 
is unsteady and episodic (see also Vorobyov & Basu 2006, 2010; Stamatellos et al. 2011). When the disk becomes 
gravitationally unstable by mass accretion or radiation cooling, non-axisymmetric structures like spiral arms are formed 
and transport mass and angular momentum, then the disk is stabilized again unless it is too unstable and fragments. 
This stochastic accretion causes the first core flickering. Its time scale is corresponding to the orbital time scales, 
ranging from a few years to a few decades. We will discuss its observational properties in a separate paper in the near 
future. 

Similar to the case with Ohmic dissipation only, the magnetic field lines in the outer region of the first core are 
strongly straightened out, running almost vertically, by strong Ohmic dissipation and ambipolar diffusion (Figures 10 
and 11). The field lines are wound up tightly near the center of the first core where the magnetic fields are coupled with 
the gas, but again the magnetic pressure is not dominant during the first core phase and the fields behave passively in 
this phase. 

Although we stop the simulation soon after the second collapse begins (when T c ~ 2,100K), we can see that the 
evolution during the second collapse is already different from the other two models. Because it is supported by the 
centrifugal force, the gas does not collapse dynamically even when it loses the gas pressure support. It rather collapses 
gradually with a time scale much longer than the free fall. Unfortunately, we have to stop the simulation at this point 
because the Jeans condition is violated and fragmentation occurred. We cannot tell but only can speculate how the 
system will evolve later, but it is likely that a compact binary or multiple system is formed (Machida et al. 2008b). 

To summarize, ambipolar diffusion is so effective and drastically changes the first core structure. It removes magnetic 
flux from the first core, suppresses angular momentum transport, and enables early formation of rotationally supported 
disks, even before formation of the protostar. The disk is massive enough to become gravitationally unstable and 
becomes non-axisymmetric. On the other hand, the disk size remains small in the early phase, contrary to purely 
hydrodynamic models without magnetic fields which yields very large disks supported by rotation even in the first core 
phase (Bate 1998, 2010; Saigo et al. 2008; Tsukamoto & Machida 2011; Tomida 2014). Thus, these non-ideal MHD 
effects can solve the magnetic braking catastrophe and maintain consistency with observations of young stellar objects 
(e.g. Maury et al. 2010). 


5. DISCUSSION 

5.1. Comparison with the Precedent Non-Ideal MHD Simulations 

Our results seem to be inconsistent with the earlier non-ideal MHD simulations by Li et al. (2011) (see also Mellon & 
Li 2009), in which no rotationally supported disk was formed even with ambipolar diffusion. Because there are many 
differences in the simulation setups including the initial conditions, comparison is not trivial but at least the diffusion 
rates are similar. We suspect that the most important difference between their simulations and ours is the boundary 
condition at the center of the cloud. They have a relatively large hole with a radius of 6.7 AU, while we do not have 
any boundary there but resolve the small-scale structure directly using the nested-grid technique. In fact, because the 
disks formed in our simulations are small, it is not inconsistent with their simulations, at least in the early phase of 
protostellar collapse. 

Our calculations are expensive and we cannot follow the long-term evolution, but this high resolution is required to 
resolve the small first core disk. We believe it is crucial to resolve the first core because the non-ideal MHD effects 
work effectively and a rotationally supported disk is formed there. Otherwise, these non-ideal MHD effects have only 
minor influence on the large scale evolution. Indeed, Machida et al. (2014) demonstrated that the sink particle or inner 
hole has a significant impact on the results. We propose that the effects of these numerical tricks can be estimated by 
measuring the angular momentum accreted into the inner boundary or sink particles. 

5.2. Perspectives of Future Observations 

As discussed in the introduction, recent interferometric observations of Class-0 objects, the youngest protostars 
including L1527 IRS, suggest that circumstellar disks in the early phase of star formation should be small, i.e. 
.RlSlOOAU, while they can be formed in the early phase (Maury et al. 2010; Ohashi et al. 2014). These obser¬ 
vational results do not favor previous hydrodynamic simulations of protostellar collapse without magnetic fields which 
tend to produce large rotationally supported disks in the early phase (e.g. Bate 1998, 2010; Saigo et al. 2008; Tsukamoto 
& Machida 2011). In contrast, our model with both ambipolar diffusion and Ohmic dissipation is consistent with this 
picture. In order to test the scenario, more observations of very young stellar objects with higher resolutions that can 
sufficiently resolve the circumstellar disks are required. In particular, if we could directly observe a first core or at least 
its remnant shortly after the second collapse, it would be a critical test for our models. If we could detect a compact 
rotating disk in a collapsing cloud but associated with no significant stellar signature (e.g. near infrared emission or 
high velocity jets), it would be crucial evidence indicating that a rotationally supported disk is already formed in the 
first core phase. 

From numerical simulations (e.g. Boss & Yorke 1995; Tomida et al. 2010a; Tomisaka & Tomida 2011; Saigo & 
Tomisaka 2011; Commergon et al. 2012), a first core is expected to (1) have a spectral energy distribution like a low- 
temperature gray body without a hot infrared stellar emission, (2) be a compact, condensed core which may exhibit 
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log (Gas Density [g cnr 3 ]) 

Fig. 10.— Magnetic field lines (yellow lines) threading the first cores and surrounding pseudo-disks (1=11 or ~ 18 AU) at the end of the 
first core phase. Vertical and horizontal density cross sections are also presented. Note that the density of the field lines does not reflect 
the strength of the magnetic fields. 
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Fig. 11.— Vertical cross sections of the magnetic Reynolds numbers in the first cores and surrounding envelopes (1=10 or ~ 35AU), when 
the temperature at the center of the cloud is 600K (left column), 1,400K (middle) and 2,000K (right). We use the gas velocity and local 
Jeans length to evaluate the Reynolds numbers. The top row shows the Reynolds number calculated from the Ohmic resistivity in Model 
O. The bottom two rows are Model OA, while the middle row indicates the Ohmic Reynolds number and the bottom row is calculated 
using the ambipolar diffusion rate. The non-ideal MHD effects work effectively in the blue regions while the red regions are close to ideal. 


a signature of infall motion, and (3) be associated with slow compact molecular outflows, but no fast jets. Although 
first cores should be very rare due to the short lifetime that only one out of a few hundreds of molecular cloud cores 
are expected, some candidates are already reported (Chen et al. 2010; Enoch et al. 2010; Chen et al. 2012; Pezzuto 
et al. 2012; Hirano & Liu 2014; Tokuda et al. 2014). Among them, L1451-nnn (Pineda et al. 2011) is one of the best 
examples which satisfies all the above criteria. We expect that direct observation of a first core or its remnant with 
high resolution interferometers, especially ALMA, will be achieved in the near future. 

5.3. Non-Ideal MHD Effects 

While the non-ideal MHD effects can avert the magnetic braking catastrophe, the magnetic flux problem may remain 
a serious problem. As shown in Figure 6, the non-ideal MHD effects increase the mass-to-flux ratio of the first core but 
only by one order of magnitude. This is insufficient to explain the difference of magnetization between molecular cloud 
cores and protostars. This is partly because we do not consider shielding of cosmic rays (see also Dapp et al. 2012; 
Padovani et al. 2014), while the ionization rate in the dense region would largely vary depending on the abundance 
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of radioactive nuclei as discussed in 2.3.2. Possibly, higher non-ideal MHD rates or other mechanisms to accelerate 
diffusion (e.g. turbulence) are required. If we use higher rates, the disk size will be also larger. Or, we speculate that 
further redistribution of magnetic fields from the formed protostar will take place in the later phase (e.g. Braithwaite 
2012 ). 

The non-ideal MHD effects play an important role in numerical simulations; they actually help convergence. In the 
ideal MHD approximation, dissipation of magnetic fields is totally determined by resolution. In protostellar collapse 
simulations, the mid-plane is a current sheet where the magnetic fields above and below are anti-parallel. In ideal 
MHD, the current sheet can be extremely thin and convergence is not well defined, at least difficult to achieve. For 
example, Joos et al. (2012, 2013) reported poor convergence of their ideal MHD simulations with nhs-aligned magnetic 
fields or turbulence; more than 20 cells per the local Jeans length are required for the nhs-aligned cases while 15 cells 
seem to be sufficient for the aligned cases. With the non-ideal MHD effects, on the other hand, the current sheet has a 
finite thickness, and therefore convergence is now well defined and can be achieved with a finite resolution. However, 
it still requires high resolution to resolve the dissipation scales sufficiently because the dissipation scale is very small 
due to the small diffusivity in the low density region. 

Although our simulations include as many physical processes as possible, a possibly important effect missing in this 
work is the Hall effect. The Hall effect is another non-ideal MHD effect which can spin up the gas and help disk 
formation (Krasnopolsky et al. 2011; Li et al. 2011; Braiding & Wardle 2012). Because the rate of the Hall effect is as 
high as the ambipolar diffusion (Wardle & Ng 1999; Nakano et al. 2002; Li et al. 2011), it is likely important in star 
and disk formation processes. However, it is not easy to implement in high-resolution 3D simulations, especially when 
we resolve the first core directly where the Hall effect works most actively. This is because it requires very small time 
steps proportional to Ax 2 /r]n where /yn is the rate of the Hall effect. The STS acceleration is not effective for the Hall 
effect because of its hyperbolic nature. Moreover, unlike Ohmic dissipation and ambipolar diffusion, the Hall effect 
changes the characteristic speeds of the MHD waves (e.g. whistler waves). Therefore, to capture this effect accurately, 
it is inadequate in principle to use the operator-splitting and sub-cycling techniques for the Hall effect, like we did 
for Ohmic dissipation and ambipolar diffusion. Instead, we have to solve the MHD part and the non-ideal MHD part 
using the same, very small time steps required by the Hall effect, which will be extremely expensive. 

Finally, it should be noted that there are significant uncertainties in the rates of the non-ideal MHD effects because 
of large uncertainties in the underlying dust properties as discussed in 2.3.2. The uncertainties also exist in the dust 
opacities used in the radiation transfer part. Because these uncertainties have direct impacts on formation and evolution 
of first cores, the results, especially the disk size, might be easily changed by an order of magnitude. Therefore, this 
work should be considered as a case study based on simple assumptions on the dust properties. A broad parameter 
survey is required to clarify the effects of the dust properties such as the structure and size distribution. It will be 
important and interesting to solve evolution of dust grains, including both growth and destruction, consistently in 
parallel with RMHD simulations. There are many topics on dust left to be explored for the future works. 

6 . CONCLUSIONS 

In order to study the effects of Ohmic dissipation and ambipolar diffusion on protostellar collapse and disk formation, 
we performed 3D nested-grid RMHD simulations. We simulated only the early phase of star formation before formation 
of protostellar cores, but we found that significant differences already arose at this point. We revealed that the non¬ 
ideal MHD effects qualitatively change the structure and evolution of the first cores. Especially, ambipolar diffusion 
expands the diffusive region and strongly suppresses angular momentum transport, resulting in early formation of 
rotationally supported disks. The conclusions are summarized as follows: 

1 . Ohmic dissipation works in the high density region and suppresses angular momentum transport. As a result, a 
very small rotationally supported disk is formed after the second collapse, essentially in parallel with formation 
of the protostellar core. 

2. Ambipolar diffusion extends the diffusive region to the lower density gas, and suppresses angular momentum 
transport more effectively. The first core is already supported by rotation significantly from the beginning, 
although the thermal pressure is still important. 

3. The formed disks remain small in the early phase of star formation and the rotation profile is not Keplerian yet, 
but they will grow later as gas accretion with larger angular momenta continues. 

4. The first core disk with ambipolar diffusion acquires sufficiently large mass and angular momentum to be¬ 
come gravitationally unstable and form non-axisymmetric structures. As a result, the accretion becomes nat¬ 
urally episodic. Although this is speculative at this point, it is likely that disk fragmentation occurs and a 
binary/multiple system will be formed later. 

5. Because of additional support by rotation, the first core is more massive, longer-lived, and more luminous in the 
model with ambipolar diffusion. This increases the probability of first core detections, although it is still rare. 
Also, the first core can exhibit variability by the gravitational instability with the short time scale on the order 
of years. 

These results can reconcile the observational results and theoretical studies. In previous hydrodynamic simulations 
of protostellar collapse without magnetic fields, large rotationally supported disks are formed in the early phase of 
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star formation, but those results are not favored by the observations of very young stellar objects. In ideal MHD 
simulations, on the other hand, the magnetic braking is so efficient that a circumstellar disk is difficult to form, which 
is referred as the magnetic braking catastrophe. Our results suggest that the non-ideal MHD effects can resolve the 
magnetic braking catastrophe in the early phase of star formation even before a protostellar core is formed, and enable 
formation of a small disk in the early phase. It is expected that the disk will grow later as larger angular momenta 
accrete. 

In previous studies of disk evolution and planet formation, simplified models of isolated, already-established proto¬ 
planetary disks like the minimum mass solar nebula (Hayashi 1981) have been often used as the initial and boundary 
conditions. However, our results suggest that a seed of a circum“stellar” disk can be formed in the early phase even 
before a protostar is formed, it co-evolves with the central protostar, and disk fragmentation can happen while the 
disk and star are still growing via accretion. This means that we need to consider more realistic situations in order 
to understand disk evolution and planet formation. Ultimately, we need a consistent scenario of star, disk and planet 
formation - or stellar system formation. 
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APPENDIX 


A. NUMERICAL METHODS FOR OHMIC DISSIPATION AND AMBIPOLAR DIFFUSION 

A.l. Conservative Form 

In this appendix, we describe how to update the induction and energy equations related to Ohmic dissipation and 
ambipolar diffusion. We adopt the operator-splitting approach and solve these effects separately from the other physical 
processes. The equations solved in this part are as follows: 


de 

9t 


dB 

~dt 

diss 


diss 

+ v 


■f V X j TjQ J + 
VA 


VA 


VoF ' |B | 2 

JeVxB. 


B x F | =0, 
= 0 , 


|B | 2 
(B x F) x B 

FeJxB. 


We can rewrite eq. (Al) into a conservative form: 


dB 

~dt 


+ V-T diss = 0, 

diss 


(Al) 

(A2) 


(A3) 


where T diss = T° + T A , T° and T A are the flux tensors of Ohmic dissipation and ambipolar diffusion, respectively. 
The fluxes of B b component in a-direction are: 

T? b = -no(d a B b -d b B a ), (A4) 

Ta b = ~^i(BaF b -B b F a ). (A5) 


Note that the dot product (•) of a vector and a tensor is defined as summation over the first suffix of the tensor; 
v • T = v a T ab using the Einstein summation convention, which yields a vector. With these fluxes, considering that 
these flux tensors are anti-symmetric and their diagonal components are zero (T 00 = 0), eq. (A2) is also rewritten 
like: 


de 

Ot 


+ V • (T diss • B) = 0. 

diss 


(A6) 


Again note that the dot product (•) of a tensor and a vector is summation over the second suffix of the tensor; 
T • v = T ab v b . 
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A.2. Correction Terms for divB 

Magnetic fields are defined at the cell centers in our scheme. In order to satisfy the solenoidal constraint (eq.5), we 
use the mixed cleaning scheme (Dedner et al. 2002) in the MHD part. By introducing additional correction terms for 
the non-ideal MHD part, we can improve the nature of the equations and reduce the divergence error. For Ohmic 
dissipation, this is achieved using the method proposed by Graves et al. (2008), adding correction terms proportional 
to divB in the diagonal components of the Ohmic dissipation fluxes: 

oV-B. (A7) 

For the ambipolar diffusion part, this is done by introducing additional source terms (Duffin & Pudritz 2008): 


<9B 

~dt 


corr 
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(A8) 

(A9) 


A.3. Discretization 

Discretization of eqs.(Al) and (A2) is straightforward: 


dt 


diss 


de i ’ j ’ k 


dt 
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Ax 
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(All) 


where i,j and k are the spatial indices of the cell centers in the x,y and z direction, and i ± 1/2 means the location 
of the cell surface. The Einstein notation is used in eq.(All); T xa B a = Yl a = x y z TxaB a , etc.. The fluxes at the cell 
surface (i + 1/2, j, k) are calculated using the following relations and their permutations: 


TDi-\-l,j,k _ TDi,j,k 
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(A13) 

(A14) 

(A15) 

(A16) 


For the source terms in eqs.(A8) and (A9), the derivatives at the cell centers are discretized using the following formula 
and its permutations: 


(d x B x y^ k - 


2Ax 


(A17) 


A.4- Nested-Grid Boundaries 

At the boundaries between different nested-grid levels, we have to connect solutions in different levels. Because our 
scheme uses the conservative form, this is straightforwardly implemented using the following procedure. 
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2. Calculate fluxes in all the grid levels. 

3. Recalculate the fluxes at the level boundaries in the coarser grid using the sum of the fluxes in the finer grid; 

e.g. ASFx +1/2,J ’ K ’ 1 = A S ^ +1 / 2 ’ : >^ 1+1 +F^ +1/2j+1 ’ k ’ l+1 + F l x +1/2J ’ k+1 ’ l+1 + K +1/2J+hk+hl+1 ) where (I, J,A') 
and ( i,j, k) are the corresponding cell indices in the coarser and finer grid levels, l is the index of the nested-grid 
level (larger is finer), and AS and As are the surface areas of the coarser and finer cells, respectively. 

4. Update the solution to t = t n+1 in all the grid levels using the corrected fluxes. 

Note that no time-interpolation is required because we use shared time stepping. 


A.5. Time Integration: Super Time Stepping 

As in Paper I, we adopt the Super Time Stepping method (STS, Alexiades et al. 1996; Choi et al. 2009; Connnergon 
et al. 2011) in order to accelerate the calculation when the time step for the non-ideal MHD part is much smaller than 
the MHD part. There are two parameters in STS; v which is a small positive parameter controlling the stability and 
efficiency of the scheme, and IVsts which is the number of the sub steps used in the scheme. We use v = 0.01 and 
IVsts = 6 as in Paper I. Acceleration by a factor of ~ 4 is achieved with these parameters compared to the first-order 
explicit scheme. 

We use the same time step in all the grid levels. The time step for the non-ideal MHD part is estimated from the 
standard CFL condition for the diffusion-like equations: 

At diss = 0.1 x min f —^^ , (A18) 

\Vo + VA J 

where 0.1 is a safety factor. When this time step is significantly smaller than that for the MHD part (Af d iss < 
AfMHD/AsTs)) we update the non-ideal MHD part using STS and sub-cycling. If this condition is not satisfied but 
the time step is still smaller than that for the MHD part (AtMHD/AsTS < Atdiss < A^hd), the standard forward 
Euler integrator with sub-cycling is used. If it is larger than the MHD time step (AfMHD < Af d j ss ), we update the 
non-ideal MHD part using the time step for the MHD part. 


A.6. Tests: Oblique C-shock Problems 

Here we demonstrate the validity of our non-ideal MHD solver based on the non-isothermal oblique C-shock tests 
(Wardle 1991; Mac Low et al. 1995; Duffin & Pudritz 2008; Masson et al. 2012). These are simple one-dimensional 
shock-tube tests with the non-ideal MHD effects. 

Although these problems are one-dinrensional, we make them two-dimensional tests by adopting the skew-periodic 
boundary condition. We rotate the shock tube tests by an angle of 6. When r = tanf? is a rational number, we can 
construct the skew-periodic boundary condition which does not introduce any artificial effect. In the ^-direction, we 
adopt the fixed boundary condition, which may violate the solenoidal constraint but the divB error is cleaned by the 
mixed cleaning scheme (Dedner et al. 2002) and does not affect the calculations. In the y-direction, we copy all the 
variables from the other side of the computational domain with a shift; 


■)(i, N y + 1) = v(i - rN y , 1), 


(A19) 


where N y is the number of cells in the y-direction. If we set N y so that rN y is an integer, this boundary condition 
produces no artificial effect and can be used to run one-dimensional shock tube tests on a two-dimensional grid. 

We compare the results between three configurations: fiducial liigh-resolution aligned models (Ax = 1/128 and 
r = 0), high-resolution tilted models (Ax = 1/128, r = 3/4 and N y = 4), and low-resolution tilted models (Aa: = 1/64 
and r = 3/4). The rotation angle in the tilted models is 9 ~ 36.87 deg. Each cell is square shaped; Ay = Ax. The 
resolution in our high-resolution aligned model is corresponding to the highest AMR level in Masson et al. (2012). 

We adopt the same parameter sets as Masson et al. (2012). For Ohmic dissipation C-shock, the initial left state is 
(p,v x ,v y ,B x ,B y ,P) = (0.4,3, 0,^/2, v / 2/2,0.4), the right is (p, v x , v y , B x , B y , P) = (0.71084,1.68814,0.4299, V^/2, 
1.43667,1.19222), and the uniform resistivity rj o =0.1 is used. For ambipolar diffusion, the left is (p, v x ,v y , B x , B y , P) 
= (0.5,5, 0, s/2, s/2, 0.125), the right is (p, v x , v y , B x , B v , P ) = (0.9880, 2.5303,1.1415, s/2, 3.4327,1.4075), and the am¬ 
bipolar diffusion rate is density dependent; = l/(75p). We run the simulations until steady solutions are obtained. 
In both calculations, the adiabatic index is T = 5/3, the CFL number is 0.25, and the STS time integrator is selected. 
The results are shown in Figure 12. Overall, we successfully obtained good numerical solutions consistent with the 
senri-analytic solutions. In particular, we found no problem in the multi-dimensional (i.e. tilted) models. 
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Fig. 12.— The results of the oblique C-shock tests with Ohmic dissipation (left) and ambipolar diffusion (right). 
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